#include "Cell.h"

const double M[Q_DIM][Q_DIM] =		{{1, 1, 1, 1, 1, 1, 1, 1, 1},
										{-4,-1,-1,-1,-1, 2, 2, 2, 2},
										{ 4,-2,-2,-2,-2, 1, 1, 1, 1},
										{ 0, 1, 0,-1, 0, 1,-1,-1, 1},
										{ 0,-2, 0, 2, 0, 1,-1,-1, 1},
										{ 0, 0, 1, 0,-1, 1, 1,-1,-1},
										{ 0, 0,-2, 0, 2, 1, 1,-1,-1},
										{ 0, 1,-1, 1,-1, 0, 0, 0, 0},
										{ 0, 0, 0, 0, 0, 1,-1, 1,-1}}; //M[i][j]

const double Minv[Q_DIM][Q_DIM] ={{ 1./9., -1./9.,  1./9.,     0,      0,     0,      0,     0,    0},
										{ 1./9.,-1./36.,-1./18., 1./6., -1./6.,     0,      0, 1./4.,    0},
										{ 1./9.,-1./36.,-1./18.,     0,      0, 1./6., -1./6.,-1./4.,    0},
										{ 1./9.,-1./36.,-1./18.,-1./6.,  1./6.,     0,      0, 1./4.,    0},
										{ 1./9.,-1./36.,-1./18.,     0,      0,-1./6.,  1./6.,-1./4.,    0},
										{ 1./9., 1./18., 1./36., 1./6., 1./12., 1./6., 1./12.,     0, 1./4.},
										{ 1./9., 1./18., 1./36.,-1./6.,-1./12., 1./6., 1./12.,     0,-1./4.},
										{ 1./9., 1./18., 1./36.,-1./6.,-1./12.,-1./6.,-1./12.,     0, 1./4.},
										{ 1./9., 1./18., 1./36., 1./6., 1./12.,-1./6.,-1./12.,     0,-1./4.}}; //M^{-1}[i][j]

Cell::Cell(void) {
	for (int l=0; l<Q_DIM; l++) {
		f[l] = 0.0;
		m[l] = 0.0;
		s[l] = 0.0;
		meq[l] = 0.0;
	}
}

void Cell::updateRhoAndU(void) {
	double _dr = 0.0;
	double _jx = 0.0;
	double _jy = 0.0;
	for(int j=0; j<Q_DIM; j++) {
		_dr += M[iDR][j]*f[j];
		_jx += M[iJX][j]*f[j];
		_jy += M[iJY][j]*f[j];
	}
	dr = _dr;
	ux = _jx/r;
	uy = _jy/r;
}

void Cell::calcMomenta() {
	for(int i=0; i<Q_DIM; i++) {
		m[i] = 0.0;
		for(int j=0; j<Q_DIM; j++)
			m[i] += M[i][j]*f[j];
	}
}

void Cell::calcDF() {
	for(int i=0; i<Q_DIM; i++) {
		f[i] = 0.0;
		for(int j=0; j<Q_DIM; j++)
			f[i] += Minv[i][j]*m[j];
	}
}

void Cell::calcEQ() {
	meq[iE] = -2.0*alpha2*m[iDR] + 3.0*(m[iJX]*m[iJX] + m[iJY]*m[iJY]);
	meq[iEPSILON] = m[iDR] - 3.0*(m[iJX]*m[iJX] + m[iJY]*m[iJY]);
	meq[iQX] = -m[iJX];
	meq[iQY] = -m[iJY];
	meq[iPXX] = m[iJX]*m[iJX] - m[iJY]*m[iJY];
	meq[iPXY] = m[iJX]*m[iJY];
}

void Cell::collide() {
	calcMomenta();
	calcEQ();
	for (int j=1; j<Q_DIM; j++)
		m[j] = m[j] - s[j]*(m[j] - meq[j]);
	calcDF();
	updateRhoAndU();
}

void Cell::collide(double _ux, double _uy) {
	calcMomenta();

	meq[iJX] = _ux*r;
	meq[iJY] = _uy*r;
	m[iJX] = _ux*r;
	m[iJY] = _uy*r;

	calcEQ();

	for (int j=1; j<Q_DIM; j++)
		if ((j != iJX) && (j != iJY))
			m[j] = m[j] - s[j]*(m[j] - meq[j]);

	calcDF();
	updateRhoAndU();
}

void Cell::setU(double _ux, double _uy, double _r) {
	m[iDR] = 0.0;
	m[iJX] = _ux*_r;
	m[iJY] = _uy*_r;
	r = _r;

	calcEQ();
	m[iE] = meq[iE];
	m[iEPSILON] = meq[iEPSILON];
	m[iQX] = meq[iQX];
	m[iQY] = meq[iQY];
	m[iPXX] = meq[iPXX];
	m[iPXY] = meq[iPXY];

	calcDF();
	updateRhoAndU();
}
